Valence-band mixing in first-principles envelope-function theory 
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This paper presents a numerical implementation of a first-principles envelope-function theory de- 
rived recently by the author [B. A. Foreman, Phys. Rev. B 72, 165345 (2005)]. The examples studied 
deal with the valence subband structure of GaAs/AlAs, GaAs/Alo.2Gao.sAs, and Ino.53Gao.47As/InP 
(001) superlattices calculated using the local density approximation to density-functional theory and 
norm-conserving pseudopotentials without spin-orbit coupling. The heterostructure Hamiltonian is 
approximated using quadratic response theory, with the heterostructure treated as a perturbation of 
a bulk reference crystal. The valence subband structure is reproduced accurately over a wide energy 
range by a multiband envelope-function Hamiltonian with linear renormalization of the momentum 
and mass parameters. Good results are also obtained over a more limited energy range from a 
single-band model with quadratic renormalization. The effective kinetic-energy operator ordering 
derived here is more complicated than in many previous studies, consisting in general of a linear 
combination of all possible operator orderings. In some cases the valence-band Rashba coupling 
differs significantly from the bulk magnetic Luttinger parameter. The splitting of the quasidegen- 
erate ground state of no-common-atom superlattices has non-negligible contributions from both 
short-range interface mixing and long-range dipole terms in the quadratic density response. 

PACS numbers: 73.21.-b, 73.61.Ey, 71.15.Ap 
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I. INTRODUCTION 

In a recent paper, the author has developed a first- 
principles multiband envelope-function theory for semi- 
conductor heterostructures i This theory can in princi- 
ple provide accurate numerical predictions of envelope- 
function Hamiltonian parameters, but only if a reli- 
able quasiparticle self-energy is used as input. On the 
other hand, if the input potential is taken from a sim- 
ple density-functional calculation^ the numerical values 
are less accurate, but the theory can still provide deep 
insight into the basic physics of the interface and clarify 
various limitations of commonly used envelope-function 
models. The purpose of the present article is to ex- 
plore these topics using numerical examples taken from 
density-functional calculations of the valence subband 
structure in semiconductor superlattices. 

In order to obtain the most direct comparison with 
envelope-function theory as it is commonly used, the 
Hamiltonian was derived in Ref. [l| in the form of a matrix 
containing differential operators and energy-independent 
functions of position. Using a differential equation of low 
order to describe an abrupt heterojunction might seem 
at first glance to be a gross theoretical blunder, since it is 
well known that bulk effective-mass theory 3 - is valid only 
for slowly varying perturbing potentials. Yet effective- 
mass theory is widely accepted as a valid lowest-order ap- 
proximation for the shallow impurity problem^££ even 
though the atomic impurity potential is not slowly vary- 
ing. Its validity in this context is established by the ex- 
istence of rigorous methods for extending effective-mass 
theory to include higher-order perturbations^^ includ- 
ing the rapidly varying part of the impurity potential. 8 
These produce short-range and long-range corrections to 
the effective Hamiltonian (each with its own independent 
parameters), which lead to the well known chemical shift 



and splitting of effective-mass degeneracies^^ thereby 
bringing the rough estimates of effective-mass theory into 
much closer agreement with experiment. 

A heterostructure is nothing but an assembly of atoms. 
It is therefore difficult to imagine a valid argument for 
accepting the use of generalized effective-mass differen- 
tial equations in the shallow impurity problem, yet cat- 
egorically rejecting them for small atomic perturbations 
in heterostructures. According to linear response the- 
ory, the perturbation generated by a collection of small 
atomic perturbations is to leading order just the super- 
position of the individual perturbations. The rapidly 
varying part of the heterostructure potential 9 -^ can thus 
be handled by the same techniques used in the impurity 
problem^ 8 - The validity of low-order differential equations 
for shallow heterostructures — as for shallow impurities — 
consequently rests on two fundamental assumptions: 3,8 
that the envelope functions satisfying these equations are 
slowly varying 9 (i.e., have a Fourier transform limited to 
small wave vectors) within the energy range of experi- 
mental interest, and that the atomic perturbations are 
truly "shallow" in real heterostructures. 

The first assumption is confirmed by numerical exam- 
ples in later sections of this paper, in agreement with 
prior empirical pseudopotential studiesi 11 ' 12 ! 13 ! 14 For the 
second assumption, the definition of "shallow" is a rela- 
tive one that depends on the energy separation between 
the states included explicitly in the envelope-function 
model and the remote states treated as perturbations 
[see, e.g., Eq. (11.33) of Ref. 0|. Consider the case of a 
GaAs/AlAs heterostructure, where the valence-band off- 
set is about 0.5 eV and the energy gap is about 1.5 eV. 
Treating GaAs/AlAs as a shallow perturbation of the 
virtual crystal Alo.5Gao.5As would be expected to yield 
marginal results (at best) in a single-band effective-mass 
model for the degenerate T valence states. Such a model 
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would be expected to give good results only for weaker 
perturbations, such as those in a GaAs/Al 2 Gao.sAs het- 
erostructure (within the virtual-crystal approximation). 
On the other hand, treating GaAs/AlAs as a shallow 
perturbation would be expected to produce good results 
in a multiband envelope-function model that includes a 
few of the nearest conduction-band states explicitly, since 
the remote bands are then more than 5 eV away from the 
valence-band maximum. (This method of making moder- 
ately "deep" perturbations "shallow" by including more 
bands in the Hamiltonian was proposed by Keldysh for 
the deep impurity problem in bulk semiconductors. 15 ) All 
of the above expectations are confirmed below by numer- 
ical examples for GaAs/Alo.aGao.gAs, GaAs/AlAs, and 
Ino.53Gao.47As/InP. 

To be more specific, the multiband envelope-function 
Hamiltonian of Ref. 0] was derived by treating the hetero- 
structure as a perturbation (within the pseudopotential 
approximation) of a virtual bulk reference crystal. The 
self-consistent potential energy of the heterostructure 
was approximated using the linear and quadratic terms of 
nonlinear response theory (see the preceding paper— for 
further details). A finite-order envelope- function Hamil- 
tonian was constructed by using Luttinger-Kohn pertur- 
bation theory 3 - ^ 16 i 17 to eliminate the k • p and potential- 
energy coupling to remote bands, working to order fc 4 in 
the bulk reference Hamiltonian) 1 ^ to order k 2 in the lin- 
ear response termS ) 16 i 17 and to order k° in the quadratic 
response^ This theory is shown here to work well in 
a 3-state model for the Ti5 V valence states (i.e., ne- 
glecting spin-orbit coupling) of a GaAs/Alo^Gao.sAs su- 
perlattice and in a 7-state {ri 5c , Ti c , ri 5v } model for 
GaAs/AlAs and Ino.53Gao.47As/InP superlattices. Ex- 
amples illustrating the success of a 4-statc {ri c ,ri 5v } 
model for Ino.53Gao.47As/InP superlattices have been 
given elsewhere^ 

To test the limits of the single-band (ri5 V ) model, this 
paper also extends the theory of Ref. [l| to include terms of 
order 9 2 k and 8 2 k 2 , where 9 denotes the heterostructure 
perturbation. This is shown to yield much better predic- 
tions of the position-dependent effective masses of both 
GaAs/AlAs and Ino.53Gao.47As/InP than the 0{9 1 k 2 ) 
theory of Ref. [J The resulting superlattice band struc- 
tures are also of reasonably good accuracy, but over a 
more limited energy range than the multiband models. 
The principal limitation of the single-band theory seems 
to be the presence of spurious solutions generated by 
0(9°k 4 ) terms in the bulk reference Hamiltonian. 

The extended 0(9 2 k 2 ) theory also yields some interest- 
ing conclusions regarding operator ordering in the effec- 
tive kinetic-energy operator T in effective-mass theory. 
Much previous work has examined various possible ways 
of choosing the exponents a, /3, and 7 in the von Roos 
parametrization 1 ^ 

T v r = \{jn a pm^ pm 1 + m? 'pm^ ] pm a ) , (1) 

where p is the momentum operator (in one dimension), 
m is the effective mass, and a -1-/3 + 7 = —1. Morrow and 



Brownstei n 20 ! 21 have argued that only exponents satis- 
fying a — 7 are physically permissible in abrupt hetero- 
structures, which would rule out seemingly reasonable 
possibilities such as^ 2 - T = j(m~ 1 p 2 + p 2 m^ 1 ). 

The present theory is not consistent with any single op- 
erator ordering of this type. Instead, the terms of order 
p 2 derived here take the form of a linear combination of 
terms containing position-dependent functions having all 
possible operator orderings with respect to p. The appar- 
ent conflict between this result and the theory of Morrow 
and Brownstei n 20 ! 21 is resolved by the fact that these are 
smooth functions of position with no discontinuity even 
at an abrupt junction. 

The structure of the paper is as follows. The quadratic- 
response theory used in the preceding article 2 - to sim- 
plify the self-consistent pseudopotential is briefly re- 
viewed in Sec. HU Section Mil describes the construc- 
tion of the envelope-function Hamiltonian, which is also 
studied from the perspective of the theory of invariants. 
In Sec. IIV1 the parameters of the Hamiltonian are calcu- 
lated and discussed for the material systems GaAs/AlAs, 
GaAs/Alo^Gao.sAs, and Ino.53Gao.47As/InP. The va- 
lence subband structures for (001) superlattices of these 
materials are calculated in Sec. |V1 which compares 
the predictions of various approximate envelope- function 
models with "exact" numerical calculations. Finally, the 
results of the paper are discussed and summarized in Sec. 

eh 



II. QUADRATIC-RESPONSE THEORY 

The foundations for this work were laid in the pre- 
ceding article? 2 in which methods are developed for ap- 
proximating the self-consistent heterostructure pseudo- 
potential using quadratic-response theory. This section 
presents a brief summary of the main ideas. 

All of the numerical results in this paper are taken from 
self-consistent total-energy calculation s 23 ! 24 performed 
in a plane-wave basis using the ABINIT softwar o 25 ' 26 ! 27 
with nonlocal norm-conserving pseudopotentials and the 
local-density approximation (LDA) to density-functional 
theory. Spin-orbit coupling is neglected. This model 
was chosen in order to permit direct comparisons of 
envelope-function calculations with "exact" numerical 
calculations in the relatively large superlattices where 
envelope-function theory is valid. Some technical ingre- 
dients of the calculations and the justification for choos- 
ing this particular model are discussed further in Ref. 2. 
As shown there, the model system predicts valence-band 
offsets in reasonably good agreement with experiment, 
but (as usual for LDA calculations 2 —) does not predict 
accurate conduction bands. 

The envelope-function Hamiltonian in this paper is 
not calculated directly from the "exact" superlattice 
pseudopotential (since that would require a separate 
self-consistent calculation for each new structure), but 
is instead obtained from the linear and quadratic re- 
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sponse to virtual-crystal perturbations of a bulk refer- 
ence crystal*^ The reference crystal is denned as the 
virtual-crystal average of the bulk constituents (e.g., 
Alo.5Gao.5As for GaAs/AlAs and Ino.765Gao.235Aso.5Po. 5 
for Ino.53Gao.47As/InP). The perturbation of the hetero- 
structure relative to the reference crystal is defined by 
the change in pseudopotential 

AVp S p(x) = £ £ « n (x - Rq), (2) 

a R 

where w ; " n (x) is the ionic pseudopotential of atom a, R Q 
is the position of atom a in unit cell R, and # R is the 
change in fractional weight of atom a in cell R of the 
heterostructure relative to the reference crystal. 

The fundamental assumption of nonlinear response 
theory is that the valence electron density n and the 
self-consistent pseudopotential can be expressed as power 
series in the perturbation # R ; e.g., 

n(x) = n (0) (x) + (x) + n {2) (x) + • • • , (3a) 
nWfxJ^'Mx), (3b) 

Q,R 

n< 2 >(x) = Y! £' 0&&AnS£,(x). (3c) 

a,R a',R' 

Here n(°)(x) is the density of the reference crystal, 
n^^(x) is the total linear response, and 7?/ 2 )(x) is the 
total quadratic response. The coefficients An R (x) and 
A n RR'( x ) give the linear and quadratic response to in- 
dividual monatomic and diatomic perturbations of the 
reference crystal. The summations in Eqs. (I3b[) and (|Bc|) 
are limited to independent values of ai 2 - 

To find all eigenstates of the Hamiltonian, one must 
know the electron density n(k + G) for all values of the 
crystal momentum k. However, the properties of slowly 
varying envelope functions depend on the density only 
in a small neighborhood of each bulk reciprocal-lattice 
vector G. Within such a neighborhood, it is convenient 
to approximate n(k+G) using a power-series expansion. 
In a superlattice, the allowed values of k + G satisfy 
k = kz for small |k|, where z is the direction normal to 
the interface plane. The power series for the monatomic 
and diatomic response coefficients in Eqs. (|3b[) and (j3"c|) 
therefore have the form 

00 

An(k + G) =J2(~ik) l Ani(G). (4) 

1=0 

Similar expansions are used for the local part of the ionic 
pseudopotential and the exchange-correlation potential. 
In the present work, the power series for the potentials 
are approximated using < / < 2 for the linear potential 
and I = for the quadratic potential^ 

The nonlocal pseudopotential is also expanded in a 
similar power series, retaining terms of order fc 4 in the 
reference crystal and order k 2 in the heterostructure per- 
turbation (which is linear by definition). The methods 



used to obtain this expansion are described in Appendix 

El 

III. ENVELOPE-FUNCTION HAMILTONIAN 

The k-space expansion coefficients for the linear and 
quadratic density and short-range potentials may now 
be used to construct an envelope- function Hamiltonian. 1 
The entire set of expansion coefficients including the local 
and nonlocal potentials is subjected to the unitary trans- 
formation given in Eq. (4.15) of Ref. l", which converts 
all matrices from the plane-wave basis to a Luttinger- 
Kohn basis of zone-center Bloch functions for the refer- 
ence crystal. Perturbation theory is then used to reduce 
these 283 x 283 matrices (where 283 is the number of 
r plane waves in the reference-crystal basis 2 -) to a much 
smaller dimension by eliminating the coupling between 
states in the set A of interest (e.g., A = {Pi5 V }) and all 
other states. The Luttinger-Kohn unitary transforma- 
tion methoe£ & 16 ' 17 i 28 is used in order to obtain energy- 
independent Hamiltonian parameters. The basic pertur- 
bation formulas are summarized in Appendix [Bl 

The renormalized Hamiltonian for set A in the ref- 
erence crystal is defined by its matrix elements (in the 
Luttinger-Kohn basis) 

(nk\H^\n'k') = {E n 5 nn , + fc,;< m , + kikjD^, 

+ hkjkiC^l, + k i k ] kik m Q' b ^J l )8\ skl , (5) 

the coefficients of which are defined in Appendix[Cj Since 
the set A = {ri5 V } will be the focus of subsequent numer- 
ical study, symmetry restrictions on the kinetic momen- 
tum matrix Tr l nn , and the inverse effective-mass tensor 
fl^, are of interest. From the symmetry of the zinc- 
blende crystal, the matrix tt 1 must have the form 

7T* = -iRlctffc |{Vfc}, (6) 

where ey^ is the antisymmetric unit tensor, {AB} = 
^(AB + BA), and Ij is the j component of the orbital an- 
gular momentum in a basis A = {\X), \ Y), \ Z)} of p-like 
orbitals<22, The coefficient R = -i7r| y is real by time- 
reversal symmetry, and the hermiticity of the Hamilto- 
nian requires R = — R = 0. Likewise, the matrix is 
given by22 

l>" = LSijl + (M - l...i,jf - N(l - <),; .{ 1, 1,\ 

+ iKe ijk I k , (7) 

where 1 represents the 3x3 unit matrix and the constants 
L = D% X ,M = £>»V N = D x x y Y + D y x x Y , and K = 
D^y- — D V XY are real. Here L, M, and N determine the 
anisotropic Ti5 V effective masses, while K is the effective 
Lande g factor for the valence band.— 

The renormalized linear-response Hamiltonian for set 
A is given by^ 

(nk\BW\n'k') = ]T'tf Q (k-k')iO(k,k'), (8) 

a 
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where 9 a (k) is the Fourier transform^ of the atomic dis- 
tribution function 9^ and 

iO(k,k') = E% n , + E a 5 nn ,6w + kA + Kn>K 

Here S Q = Ann^O) is the nonanalytic contribution^^ 
from the linear quadrupole moment of the electron den- 
sity, which merely shifts the mean energy by a constant; 
the remaining terms are defined in Appendix |Pl The no- 
tation has been modified with respect to Ref. Ill in order 
to draw attention to similarities with the bulk Hamil- 
tonian ([5]) and avoid undue proliferation of symbols. It 
should be noted that the superscripts on the coefficients 
in Eq. indicate how the coordinate and momentum 
operators are ordered; for example, in the coordinate rep- 
resentation, the term proportional to D l has the form 
D^,pi9 a (x)pj , where p is the momentum operator.— 

The symmetry of the coefficients in Eq. ([9]) is the same 
as the symmetry of site a in the reference crystal. For 
a zinc-blende reference crystal, the atomic sites have the 
same point group as the reference crystal, but in non- 
symmorphic crystals (e.g., diamond) the site symmetry 
is lower than the crystal symmetry 4 Thus, the linear mo- 
mentum matrix for A = {r 15v } has the same form as in 
bulk: 

TT ia = -iR a \e ijk \{I j I k }, (10a) 
ir a * = -iR a -\e ljk \{I 3 I k }, (10b) 

where the superscript dots are just placeholders to in- 
dicate where p is positioned with respect to Q (x). As 
in the bulk case, R' a — —i^xy ancl R a — ~^ 7T xy are 
real by time-reversal symmetry; however, for the linear 
response, hermiticity [i.e., 7r"*, = (7r^? n )*] requires only 
that R' a = —R a '. Therefore, unlike the bulk case, the 
linear R coefficients are not required to vanish. As dis- 
cussed in Ref. 0] (using a different notation) , the constant 
R' a generates a (5-function-like mixing of the \X) and \ Y) 
valence states at a (001) heterojunction^ This mixing 
is considered in greater detail below. 

The linear D tensor likewise has the same form as in 
bulk material: 

D ija _ £■■<*$.. l + { \l ■■ _ I. ;.,),,!; 

- N" a (l ~ /,/J + %K ' a e m I k , (11) 

with similar definitions for D %a i and D m K All of the 
coefficients L " Q , etc., are real by time-reversal symme- 
try. Hermiticity of the Hamiltonian requires that D°"?, = 
and D 1 ^, = which yields constraints of 

the form K " a = K a ", etc., but does not require any of 
the constants in pip to be zero. 

As discussed in the Introduction, for many heterostruc- 
tures the energy gap is not very large in comparison to 
the band offsets, which means that the linear approxima- 
tion for the momentum and mass terms used in Ref. 1 



is not very accurate in a single-band model. In an effort 
to learn more about the limits of single-band models, the 
perturbative renormalization of the momentum and mass 
terms was extended to terms quadratic in 9 a , with the 
results given in Appendix [EJ The resulting contributions 
can be written in the form 

(nk| AH ( V | n 'k') = ^ 9 a (k - k")9 p (k" - k') 

a, (3 k" 

x<f,(k,k",k'), (12) 

in which 

u nn' l K ' K ' K ) - a nn' + K ^nn' + K i 71 nn' + n nn' K i 
+ hk J D^ + k l D^k>+D^kik> 

+ kik'jD^P + KD^?k>> + D^k'Jk'j. (13) 

Here again the superscripts indicate the positioning of 
the associated operators, with (for example) the term 
proportional to n"^ having the form n"^,9 a (x)pi9p(x) 
in coordinate space. Hermiticity of the Hamiltonian gives 
constraints such as = (D 3 n ^)* . It should be noted 

that Eqs. (jX2j) and (| 13|) include only those quadratic con- 
tributions arising from perturbative renormalization; the 
other terms arising from direct multipolc expansions of 
the quadratic response have a different form and are con- 
sidered in Appendix [Fl 

For the Tisv Hamiltonian, the coefficients in Eq. (TT3"]) 
also have Td symmetry, so they are given by obvious 
generalizations of Eqs. (JTU]) and (fTTj) . The hermiticity 
constraints on the R coefficients are R a/3 ' = —R'^ a and 
= -R0 <* (which implies that R a a = 0), and one 
can also choose these coefficients to satisfy i? Q/3 ' = R@ a ' 
because 9 a and 9p commute. Likewise, the D coeffi- 
cients L, M, N, and K all satisfy constraints of the 
form K a @" = K~P a , K a '@' = K^ a , K a ^ = if 3 "", 
and K ' 3 - = K^ a \ and one can choose them to satisfy 

K c0- _ K fia- t0Q 

In a (001) heterostructure, where 9 a (x) = 9 a (z), the 
bulk Hamiltonian matrix elements of the form Lp 2 and 
Mp 2 are replaced (to order 9 2 ) by 

Lp 2 z - L^pl + L a {p 2 z 9 a + 9 a p 2 z ) + L a Pz 9 aPz 
+ L' a P{p 2 z 9 a 9p + Ocfipfi) + L al3 -p z 9 a 9 p p z 
+ L a ^(p z 9 aPz 9 + 9 p p z 9 aPz ) + L a ^9 aP 2 z 9,3, (14) 

where summation on a and (3 is implicit. This is a linear 
combination of all but one of the von Roos operators 1 - 
defined in Eq. (TTJ) . If the renormalization were extended 
to cubic order, we would find also a term 

L a 'P'~ < (9 a p z 9 l 3p z 9~ f + 9 lPz 9 pPz 9 a ). (15) 

Hence, the present derivation from perturbation theory 
supports not a single operator having the form |T]) with 
fixed exponents, but a linear combination of all possi- 
ble operator orderings. As mentioned in the Introduc- 
tion, this does not lead to any mathematical or physical 
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difficulties because 9 a is a smooth function of x. The 
question of whether any particular terms in the linear 
combination might happen to be negligible is considered 
in Sec. El 

In a similar fashion, bulk matrix elements of the form 
Np x p z are replaced in a (001) heterostructure by 

N PxPz - N^ PxPz + (TV- + 2N- a ) Px { Pzi 9 a } 
+ (N' af3 ' + APf«-/3-} + 2N- a/3 ) Px { Pz ,0 a 9p} 

+ (N a "P + NW>) P J aPz e 0> (16) 

where = |(iV a ^' + N^ '). Note that the 

final term proportional to 9 aPz 9p is not found in 
the usual symmctrization recipe for envelope-function 
Hamiltonians i 33 ' 34 

In a (001) heterostructure, the contribution of the R 
terms to the Hamiltonian is 

H R = -i2{I x I y }(R a [ Pz ,9 a ] +R ap [ Pz ,e a e l3 \ 

+ R a - 9 a [ Pz ,d p }), (17) 

which mixes the \X) and \Y) valence states because 



TABLE I: Luttinger parameters for several bulk compounds 
and their virtual-crystal averages. 



0-10 
-1 




(18) 



The function 9 a (z) is a smooth step-like function, hence 
i[ P z-,6a\ = d9 a /dz is localized at interfaces and has the 
form of a macroscopic average of the Dirac S function. 
The function 9 a \p z , 9p] is also localized, but the associ- 
ated term in (| 1T[) cannot be written as a simple derivative 
because R a = — R 13 ' 01 . Mixing of the type (jTTJ) has been 
studied by Ivchenko et alM 

The contribution of the K terms in a (001) hetero- 
structure is very similar: 

H K = ^{ Px Iy ~ P yI X )(K' a '[ Pz ,9 a ] 

+ (K- a ^ +K^^)[ Pz ,9 a e p ] 

+ (K a -P- -K^)9 a [ Pz ,9 p ]). (19) 

Here the operator ( Px I y — Py I x ) is analogous to the 
Rashba couplin g 17 i 35 i 36 ( Px a y — Py a x ) in the Tq conduc- 
tion band, so contributions of the form 1|19|) have been 
referred to as the valence-band Rashba coupling. 28 ' 37 i 38 
As discussed in Ref. [l|, this type of coupling was intro- 
duced in Ref. |3^ under an approximation that is equiva- 
lent to assuming that the Rashba coefficient is the same 
as the effective Lande factor K in bulk material. This 
has the advantage of reducing the number of unknown 
parameters, since K is known from magnetoabsorption 
experiments (see, e.g., Ref. |40|) . However, the bulk Lande 
factor to order 9 2 is given by 





GaAs 


Alo.5Gao.5As 


AlAs 


L 


-5.028 


-3.648 


-2.863 


iVl 


— 1.011 


— l.OZD 


1 1 OA 
— 1. 1ZD 


N 


-6.146 


-4.676 


-3.792 


K 


-2.086 


-0.993 


-0.514 




Ino. 53 Ga . 47 As Ino 


765Gao.235 AS0.5P0.5 


InP 


L 


-5.596 


-4.187 


-3.317 


M 


-1.385 


-1.253 


-1.119 


N 


-6.590 


-5.071 


-4.087 


K 


—2.797 


— 1.609 


—0.974 


where 








K a 


= K a +2K- a , 




(20b) 


K af3 


= K afj - + 2K- af3 + 




A (20c) 



K = K (0) + 9 a K a + 6 a $8K af3 , 



(20a) 



which shows that the Rashba coupling (TT9")) is generally 
independent of the bulk Lande factor. To linear order, 
replacing the Rashba coefficient with K would be a good 
approximation if 2\K " Q | < |iT a '|. As shown in SeclTVl 
this is true in some materials but not in others; hence, it 
cannot be presumed to hold true in general. 



IV. EFFECTIVE-MASS PARAMETERS 

In this section, the numerical values of the envelope- 
function parameters calculated for the model system are 
examined to see whether any general conclusions can be 
drawn regarding the structure of the Hamiltonian in Sec. 
IIII1 Values of the Luttinger parameters 2 - L, M, N, and 
K calculated for various bulk materials are listed in Ta- 
ble HI This table shows that the parameters are of order 1 
(in atomic units), and that their dependence on material 
composition is not linear. For example, the change (rela- 
tive to the reference crystal) in K is about —1 for GaAs 
but only | for AlAs. It should be noted that the bulk 
K values reported here include only the contributions 
from k • p renormalization, since the asymmetric terms 
in the nonlocal pseudopotential^i cannot be determined 
by polynomial fitting.— 

To obtain a measure of the accuracy of the linear and 
quadratic approximations, the change in effective-mass 
parameters for various bulk compounds relative to the 
reference crystal was calculated in these approximations 
using expressions of the form (|20j) and compared with 
the exact changes. The results of these calculations are 
given in Table [TTT The top part of the table considers the 
changes in GaAs and Al0.2Ga0.sAs relative to the refer- 
ence crystal Alo.1Gao.9As. Here the linear approximation 
is accurate to better than 10%, while the quadratic ap- 
proximation is accurate to better than 1%. Since the lin- 
ear changes are already a small perturbation in this case, 
the linear approximation for GaAs/Alo^Gao.sAs should 
be adequate for most purposes. 
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TABLE II: Difference between Luttinger parameters of bulk crystals and average reference crystal: Comparison of linear and 
quadratic approximations with exact differences. Error columns give per cent relative error. 









GaAs 










Al0.2Ga0.gAs 








Linear 




Quadratic 


Exact 


Linear 


Quadratic 


Exact 




Value 


Error 


Value 


Error 


Value 


Value 


Error 


Value 


Error 


Value 


AL 


-0.329710 


-7.22% 


-0.353495 


-0.526% 


-0.355363 


0.329710 


7.21% 


0.305924 


-0.526% 


0.307542 


AM 


-0.035986 


0.99% 


-0.035632 


-0.003% 


-0.035634 


0.035986 


-0.98% 


0.036339 


-0.004% 


0.036341 


AN 


-0.347108 


-6.83% 


-0.370680 


-0.501% 


-0.372548 


0.347108 


6.75% 


0.323536 


-0.498% 


0.325154 


AK 


-0.273357 


-8.70% 


-0.297525 


-0.623% 


-0.299389 


0.273357 


8.99% 


0.249189 


-0.644% 


0.250803 








GaAs 










AlAs 






AL 


-0.9968 


-27.8% 


-1.2706 


-7.93% 


-1.3800 


0.9968 


27.0% 


0.7230 


-7.85% 


0.7846 


AM 


-0.1933 


4.5% 


-0.1855 


0.24% 


—0.1850 


0.1933 


-3.6% 


0.2011 


0.27% 


0.2005 


AN 


-1.0919 


25.7% 


-1.3609 


-7.43% 


-1.4701 


1.0919 


23.5% 


0.8229 


-6.93% 


0.8841 


AK 


-0.7004 


35.9% 


-0.9828 


-10.05% 


-1.0926 


0.7004 


45.9% 


0.4179 


-12.92% 


0.4799 








Ino.53Gao.47 As 










InP 






AL 


-1.1067 


21.5% 


-1.3650 


-3.12% 


-1.4091 


1.1067 


27.2% 


0.8483 


-2.53% 


0.8703 


AM 


-0.1338 


2.1% 


-0.1319 


0.63% 


-0.1311 


0.1338 


-0.6% 


0.1358 


0.77% 


0.1347 


AN 


-1.2192 


19.7% 


-1.4758 


-2.84% 


-1.5189 


1.2192 


24.0% 


0.9626 


-2.13% 


0.9835 


AK 


-0.8778 


26.1% 


-1.1431 


-3.77% 


-1.1879 


0.8778 


38.1% 


0.6125 


-3.63% 


0.6355 



Note that M, which determines the mass of heavy holes 
along the (100) directions, is much more accurate than L, 
N, and K. This happens because the coupling of ri 5v to 
Tic does not contribute to M, 43 so the remote bands in 
M are more remote and the heterostructure perturbation 
for M is "shallower" than for the other parameters. As 
discussed in the Introduction, one can achieve a similar 
effect for L, N, and K by including Ti c in the set A~ 

For the case of GaAs/AlAs (with Alo.5Gao.5As as the 
reference crystal), the linear approximation is off by 
nearly 50% in some cases, while the quadratic approx- 
imation is accurate to within 13% for K and to within 
8% for the other parameters. The quadratic approxi- 
mation is more accurate for Ino.53Gao.47As/InP, with a 
maximum error of less than 4%. Although these results 
are not perfect, an accuracy of around 10% in a small per- 
turbation is good enough in many cases, and as shown 
in Sec. [V] the quadratic approximation for A — {Tisv} 
yields reasonably good band structures over a limited 
range of energies (although not as good as for a multi- 
band Hamiltonian) . 

Calculated values of various parameters in the ri5 V 
linear-response Hamiltonian ([9]) are listed in Table IIIII 
The linear changes in the bulk Luttinger parameters are 
determined by constants L a , M a , N a , and K a of the 
form (|20b|) . As noted below Eq. (fT9|) . the linear contri- 
bution to the valence-band Rashba coupling is just K' a ' . 

Many envelope-function calculations in the literature 
use the BenDaniel-Duke operator ordering^ i 33 i 34 i 44 i 45 in 
which it is assumed that |i' a '| > 2|i" Q |, |M Q | > 
2\M" a \, and \N-°-\ > 2\N" a \. Inspection of Table |TTT] 
shows that this is perhaps a tolerable approximation 
in some cases (e.g., light holes in GaAs/ AlAs), but it 
is a bad approximation in others (e.g., heavy holes in 
GaAs/ AlAs). It should be noted that Bastard^ and 
Burt^ both derive the BenDaniel-Duke ordering using 
variations of Lowdin perturbation theory^ which yields 



TABLE III: Linear parameters in the Tisv Hamiltonian. Here 
RC stands for reference crystal, and the labels light and heavy 
holes refer to the bulk properties in the (100) directions. 

RC Alo.5Gao.5As In .765Ga .235Aso.5Po.5 



a 




Ga 


As 


Ga 


Light hole 


L a 


-1.984 


-1.806 


-0.847 




L a - 


-1.341 


-0.586 


-1.002 




L a 


-0.321 


-0.610 


+0.077 


Heavy hole 


M a 


-0.387 


-0.329 


+0.130 




M a 


-0.039 


-0.109 


+0.093 




M" a 


-0.174 


-0.110 


+0.018 


k 2 mixing 


N a 


-2.181 


-2.074 


-0.771 




N - a . 


-1.542 


-0.550 


-1.085 




N - a 


-0.320 


-0.762 


+0.157 


Lande 


K a 


-1.399 


-1.287 


-0.993 


Rashba 


K a 


-1.372 


-0.464 


-1.089 






-0.013 


-0.411 


+0.048 


5 mixing 


R a 


-0.028 


-0.017 


-0.038 



energy-dependent mass parameters. This type of per- 
turbation theory cannot be used to draw conclusions 
about operator ordering in a Hamiltonian with energy- 
independent parameters, since Luttinger-Kohn pertur- 
bation theory is qualitatively different. A detailed com- 
parison of the two theories is outside the scope of this 
paper, but it will be noted that a direct application of 
Lowdin perturbation theory (using a power series expan- 
sion to treat the energy dependence of the denominators) 
to the present first-principles calculations yields values of 
L " Q , M" a , and N" a that are smaller than those in Table 
IIIII but still not generally negligible. 

As mentioned in Sec. IIIII it is also common prac- 
tice to estimate the Rashba coupling K' a ' by the 
approximation 3 ^ K' a ' ~ K a , which amounts to an ex- 
tension of the BenDaniel-Duke hypothesis to the anti- 
symmetric terms in the Luttinger Hamiltonian. 29 Table 
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TABLE IV: Quadratic parameters in the Tisv Hamiltonian. 
Here A" a/3 ' = 7V' Q ' 3 ' + A {q -' 3 '>, 7V Q "' 3 = A Q "' 3 + Art"'/ 3 '}, 
N { a -0-} = i( N a-e- + N - a -), and N^i 3 ^ = N^ 13 - - N '"' , 
with similar definitions for K. 



RC 


Alo.5Gao.5As 


Ino.765 


Gao.235 Aso.5 


P0.5 


(a, 13) 


(Ga, Ga) 


(As, As) 


(Ga, Ga) 


(As, Ga) 


L «li 


-1.087 


-0.633 


-0.668 


-0.262 


L . a p. 


-0.718 


-0.640 


-1.118 


-0.111 


L - a g 


+0.051 


+0.164 


+0.190 


+0.026 


L c.f> 


-0.007 


-0.081 


-0.002 


+0.010 




-0.232 


-0.120 


+0.036 


-0.235 


L . a .p 


-0.232 


-0.120 


+0.036 


+0.020 


M alS 


+0.0310 


+0.0014 


-0.0265 


+0.0129 




-0.0040 


-0.0119 


-0.1111 


+0.0009 


M ~<*8 


+0.0301 


+0.0209 


+0.0463 


+0.0065 


M 01 " 3 


-0.0004 


-0.0077 


+0.0018 


+0.0050 


M a - p - 


-0.0124 


-0.0104 


-0.0050 


-0.0066 


M^ 


-0.0124 


-0.0104 


-0.0050 


+0.0007 


N a0 


-1.074 


-0.638 


-0.679 


-0.252 




-0.965 


-0.787 


-1.195 


-0.225 


N ~a/3 


+0.069 


+0.183 


+0.243 


+0.035 




—0.248 


—0.216 


+0.031 


-0.097 




0.0 


0.0 


0.0 


-0.271 


K " 3 


-1.128 


-0.660 


-0.629 


-0.278 


K .af>. 


-0.933 


-0.746 


-0.966 


-0.217 


K - a g 


+0.009 


+0.133 


+0.151 


+0.019 




-0.213 


-0.180 


+0.035 


-0.099 




0.0 


0.0 


0.0 


-0.248 


Ra 


-0.0076 


+0.0034 


-0.0185 


-0.0034 


Ra -8 


0.0 


0.0 


0.0 


-0.0048 



IIIII shows that this is a good approximation for cationic 
perturbations in GaAs/AlAs and Ino.53Gao.47As/InP, 
but not for anionic perturbations in In 53 Gao.47As/InP. 
Hence, as noted in Ref. [l|, this approximation can only 
be relied on in general to produce a rough order-of- 
magnitude estimate of the Rashba parameter K' a ' . 

Numerical values for the quadratic renormalization 
terms in Eq. ((13J) arc listed in Tabic |lVl The bulk val- 
ues in this table are defined by expressions of the form 
(I20cj) . It should be noted that the present calculations 
on (001) supercells do not provide separate values for the 
constants N^' 13 '1, iV' Q/3 ', and N a since these terms al- 
ways appear together in the sums N' a>3 ' + N^ a '^ > and 
N a-(3 + N { a -{3-} ( the game ig true for K y Table[TV]is not 

discussed here beyond a brief comment that, although 
BenDaniel-Duke operator ordering is not a very good 
approximation in any case, it is typically better for light 
holes than for heavy holes, and better for cation pertur- 
bations than for anion perturbations. (Of course, since 
the position-dependent corrections to the bulk value of 
M are rather small, heavy-hole calculations are also less 
sensitive to the choice of operator ordering.) 

Finally, it should be emphasized that the numbers re- 
ported here are not expected to be accurate. They are 
merely intended to provide a crude picture of some of 
the qualitative features that would be found in a more 
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FIG. 1: (Color online) Energy band structure of bulk 
Alo.sGao.s As: comparison of exact calculation with 0(k 2 ) and 
0(k 4 ) k-p models, (a) 283-state k-p model; (b) 7-state k-p 
model; (c) 3-state k ■ p model. 



accurate quasiparticle calculation. 



V. SUPERLATTICE VALENCE SUBBAND 
STRUCTURE 

In this section, "exact" model calculations of the va- 
lence subband structure of (001) superlattices are used to 
evaluate the accuracy of various approximate envelope- 
function models. As a starting point, the bulk band 
structure of Alo.5Gao.5As (used as a reference crystal for 
GaAs/AlAs) is considered in Fig. [T] Part (a) shows the 
results when all 283 states of the plane-wave basis are 
retained in the k ■ p Hamiltonian. Here it makes little 
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difference whether the polynomials in the k • p Hamilto- 
nian are terminated at order k 2 or k ; both cases provide 
a good description throughout the Brillouin zone. In part 
(b), the set A = {ri 5c , Tic, ri 5v } contains 7 states (or 14 
with spi n 47 ' 48 ). The description of the band structure is 
still accurate near T, although spurious solutions within 
the energy gap do occur for both the 0(k 2 ) and 0(k 4 ) 
k ■ p models. 

Part (c) gives the results for A — {ri 5v }. Here there 
are no spurious solutions in the 0(k 2 ) effective-mass 
model, but the spurious solutions for the (9(fc 4 ) model 
occur at rather small wave vectors. The critical point at 
which the light-hole band has zero slope is about ^ of 
the distance to the X point. To prevent problems with 
spurious solutions, the envelope functions in a (001) su- 
perlattice should contain no Fourier components beyond 
this points Therefore, a superlattice with a period of 48 
monolayers^ was chosen as the standard test case, since 
this permits the inclusion of 9 envelope-function plane 
waves within the region \k z \ < -^{A-k/cl) = ^(2n/a). Pre- 
vious calculations on empirical pscudopotential models 
show that this is sufficient to achieve reasonably accu- 
rate resultsJi 1 ^ 1 ^ 

The features shown in Fig. [IJc) have a direct analog 
in the Dirac equation for relativistic electrons, for which 
the dispersion relation is 



E = y/p 2 C 2 + TO 2 C 4 



P 



P 



2m 8m 3 c 2 



(21a) 
(21b) 



Here the power series converges only when \p\ < mc. If 
the series is terminated at order p 4 , the slope of E(p) 
changes sign at p = \/2mc, which lies outside the re- 
gion of validity of the power series. Although the con- 
vergence radius of Luttinger-Kohn perturbation theory 
for the k • p Hamiltonian is not known, one would ex- 
pect it also to be of the same order of magnitude as the 
critical point where E{k) changes sign. Therefore, it is 
reasonable to choose this as the cutoff for plane-wave 
expansions^ even though it may lie slightly outside the 
convergence radius of the perturbation power series. 

The GaAs/Alo.2Gao.8As material system historically 
provided one of the first direct comparisons between ex- 
periment and effective-mass theory in heterostructures; 5 - 
Figure [5] shows the top 12 valence subbands in a (001) 
(GaAs)24(Al .2Ga .8As)24 superlattice. The exact model 
calculations are compared here with a 3-state Ti5 V 
envelope- function model based on the theory of Ref. U, in 
which terms of order k 3 and fc 4 are included only for the 
bulk reference crystal, the mass and momentum terms 
are linear in 9, and the potential is quadratic in 9. The 
envelope-function results are in excellent agreement with 
the exact calculations; the mean error in each of the first 
10 subbands does not exceed 0.1 meV. Note that the 
valence-band offset in this system is only 104 meV, which 
means that the heterostructure perturbation is indeed 
shallow^ even in a single-band Ti$ v model. These results 




Wave vector k 

FIG. 2: (Color online) Valence band structure of a (001) 
(GaAs)24(Alo.2Gao.8As)24 superlattice: comparison of exact 
and 3-state envelope-function (EF) calculations. The EF 
model is the same as that of Ref. CO, and the calculation in- 
cludes 9 plane waves (PW). Symmetry labels are denned in 
Refs. [5l] and the E axis corresponds to the bulk A axis, 
and the border of the figure in the E direction is 6% of the 
distance to the bulk X point (the border in the A or bulk E 
direction is the same physical distance). 



confirm that the theory of Ref. [H works very well in the 
weak-perturbation limit under which it was derived. 

The following examples provide a more detailed study 
of the effects of various approximations in systems con- 
taining stronger perturbations. Figure [3] shows the top 
12 valence subbands in a (001) (GaAs)24(AlAs)24 super- 
lattice. The exact model calculations are compared here 
with a 283-state envelope-function model that includes 
all zone-center Bloch functions explicitly, corresponding 
to Fig. QJa) in the bulk case. Figure [3ja) includes 25 
plane waves in the envelope-function model. The results 
here are very accurate, with an error corresponding to 
about a 2 meV shift that is nearly the same for all sub- 
bands. (To be more precise, the error in the ground state 
is +1.6 meV, which is almost the same as the +1.5 meV 
error in the bulk valence-band edge of GaAs calculated in 
Sec. VII C of the preceding paper.—) Figure (3^b) shows 
the effect of reducing the number of plane waves from 
25 to 9. There is little change for energies close to the 
valence-band maximum, but the error in subband 12 is 
significant. Note that the energy range here is much 
wider than in Fig. [21 although it still corresponds ap- 
proximately to the valence-band offset £ 

The effect of reducing set A to the 7 states in Fi5 C , 
Tic, and Tisv is shown in Fig. [4] Part (a) is just the 
multiband theory of Ref. [H, the single-band version of 
which was used previously in Fig. [21 The results here 
are even slightly better than in Fig. El^a) , which can only 
be attributed to a fortuitous cancellation of errors. In 
part (b) the 0(k 4 9 a ) terms are dropped. The results 
are still fairly accurate near the valence-band maximum; 
however, the peculiar behavior of (what should be) the 




e r a r z 

Wave vector k 



FIG. 3: (Color online) Valence band structure of a (001) 
(GaAs)24(AlAs)24 superlattice: comparison of exact and 283- 
state EF calculations, (a) 25 EF plane waves; (b) 9 EF plane 
waves. Both EF calculations use 0(k 4 ) bulk dispersion. 



ground state near Z shows that the plane-wave cutoff in 
this case is not quite sufficient to eliminate all effects of 
the 0(k 2 ) spurious solutions in Fig. [TJb). In Fig. [4jc) 
the 0(k 2 9 1 ) terms are dropped, so that the mass and 
momentum matrices are approximated by those of the 
reference crystal. Here the error becomes significant even 
at fairly small energies, which shows the importance of 
including linear terms for multiband models. 

In Figs. [5] and the set A is reduced even further to 
only the three r 15v states. Figure [5] shows the bands on 
the same energy scale as before, while Fig. [B] shows an 
expanded view of the region near the band edge. Part 
(a) uses the linear mass model of Ref. □ (the same as m 
Fig. [2]). A close inspection of the top three subbands in 
Fig. [f^a) shows evidence of the errors in the linear mass 
approximation displayed in Table [TTJ These errors are 
corrected in part (b), which includes quadratic correc- 
tions to the mass as well as (see Appendix [EJ terms of 
order A in the potential. This yields a noticeable im- 
provement, although the quantitative failure for higher 
excitations (due in large part to the use of only 9 plane 
waves) is still present. In part (c), the O(k 4 8 ) terms are 
dropped; since spurious solutions are no longer a prob- 
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FIG. 4: (Color online) Valence band structure of a (001) 
(GaAs)24(AlAs)24 superlattice: comparison of exact and 7- 
state EF calculations. The EF models are (a) O(k 4 0°+k 2 6 1 + 
fcV); (b) 0(k 2 6 l + fcV); (c) O(k 2 8 + k°6 2 ). All EF calcu- 
lations use 25 plane waves and O(0 2 ) potentials. 



lem, the number of plane waves is increased to 25. It 
can be seen that the top three subbands are still quite 
accurate under this approximation. 

Figure [7] shows the valence subband structure of a 
(001) (In .53Gao.47As)24(InP) 2 4 superlattice. Part (a) 
gives the results obtained from the original 283-state ba- 
sis. The error here is larger than in the analogous calcu- 
lation for GaAs/AlAs in Fig.[3fa); the ground-state error 
is +5.1 meV, which is close to the error of +4.5 meV in 
the bulk valence-band edge of Ino.53Gao.47As calculated 
111 Sec. VII C of Ref. The deviation from a constant 
shift of about 5 meV is negligible for the first 10 sub- 
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e r a r z 

Wave vector k 

FIG. 5: (Color online) Valence band structure of a (001) 
(GaAs)24(AlAs)24 superlattice: comparison of exact and 3- 
state EF calculations. The EF models are (a) O(k 4 6 + 
k 2 Q x + fcV), 9 PW; (b) O(k 4 6 + k 2 8 2 + k°0*), 9 PW; (c) 
0(k 2 9 2 + k°9 4 ), 25 PW. 



bands, which covers the full range of the valence-band 
offset^ 

In Fig. E^b) the basis is reduced to 7 states using the 
linear mass renormalization of Ref. [l|. The results are 
quite close to the exact calculation, although again the 
improvement over part (a) is fortuitous. This is demon- 
strated in part (c), which includes additional terms of 
order 9 2 k 2 and 8 4 k°. Most of the bands are shifted 
slightly upward, returning nearly to the result from part 
(a). The shift is mainly due to O(9 3 k ) corrections in the 
renormalized potential [see Eq. (|E11[) ]. Hence, neglect- 
ing 0(9 3 k°) terms in perturbative renormalization [Fig. 
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FIG. 6: (Color online) Expanded view of the top eight valence 
states from Fig. [5] 



[7jb)] approximately compensates for the neglect of cubic 
response terms in the original Hamiltonian. 

Figure 2 of Ref. HH shows the results of calculations 
that are the same as Fig.JTJb), but for the 4-dimensional 
set A — {Tic, ri 5v }. The results with O(6 k 4 ) terms are 
almost identical to Fig.[7Jb). The effect of dropping the 
O(8 k 4 ) terms is somewhat more significant than in Fig. 
0Jb), however. 

Finally, Fig. [8] shows the predictions of the single-band 
Ti5 V model for Ino.53Gao.47As/InP. When 0(k 4 ) terms 
are included in the bulk Hamiltonian, the critical point 
of zero slope in the light- hole (100) dispersion is closer to 
r than it was in GaAs/AlAs, so that now only 7 plane 
waves can be included in the envelope functions if one 
wishes to avoid spurious solutions. Apart from this re- 
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FIG. 7: (Color online) Valence band structure of a (001) 
(Ino.53Gao.47As)24(InP)24 superlattice: comparison of exact 
and EF calculations, (a) 283-state EF model; (b) 7-state 



FIG. 8: (Color online) Valence band structure of a (001) 
(Ino.53Gao.47As)24(InP)24 superlattice: comparison of ex- 
act and 3-state EF calculations. The EF models are (a) 



EF model, O(k 4 9 + k 2 9 x + k°6 2 ); (c) 7-state EF model, O(k 4 6 + k 2 6 4 + k°6 2 ), 7 PW; (b) O(k 4 9 + k 2 9 2 + k°9 4 



O(k 4 9 + k 2 9 2 + k°9 4 ). All EF calculations use 25 plane waves 
and 0(k 4 ) bulk dispersion. 



striction, the envelope- function models used in parts (a), 
(b), and (c) of Fig. [5] are the same as those used in the 
corresponding parts of Figs. [5] and [5] It can be seen that 
in this case the limitations of using only 7 plane waves are 
sufficiently severe that, for the first seven subbands, one 
is better off omitting the 0(k 4 ) terms and including more 
plane waves. It should be noted that the predictions of 
the single-band Hamiltonian for real Ino.53Gao.47As/InP 
superlattices would likely be substantially worse than 
what is shown here (perhaps even qualitatively incor- 
rect), since the energy gap of Ino.53Gao.47As in the model 
system is 61% larger than the experimental value (see 



7 PW; (c) 0(k 2 6 2 + k°9 4 ), 25 PW. 



Sec. IV of the preceding paper—). 

Although it is not visible on the scale of these fig- 
ures, the double degeneracy of the f ground state in 
GaAs/AlAs is removed in Ino.53Gao.47As/InP due to the 
reduction in symmetry from to Civ- The primary 
cause of the splitting is mixing of the \X) and \Y) va- 
lence states due to the short-range interface mixing in Eq. 
(jTTJ) and the long-range interface dipole potential in Fig. 
10 of the preceding paper £ The splitting of the quaside- 
generate ground state calculated exactly and in various 
envelope- function models is presented in Table [Vj It can 
be seen that all of the envelope-function models give a 
reasonably good estimate of the splitting (which means 
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TABLE V: Splitting of the F ground-state degeneracy in a 
(001) (Ino.53Gao.47As) 2 4(InP) 2 4 superlattice. 



Model 


Diatomic dipole included? 
Yes No 


Exact 


0.639 meV 




EF (283 states) 


0.722 meV 




EF (7 states) 


0.626 meV 


0.426 meV 


EF (3 states) 


0.585 meV 





that they provide a satisfactory description of the mi- 
croscopic wave function in the interface region). How- 
ever, when the diatomic dipole terms in Fig. 10 of Rcf. 
[H (and their associated polarization of the bulk reference 
crystal) are omitted, the splitting of the ground state 
in the 7-dimensional envelope-function model is reduced 
by about one third. This shows that the practice of fit- 
ting experimental splitting data to short-range interface 
terms onl y 32 ' 53 ! 54 may give an incorrect description of 
the basic physics and overestimate the magnitude of the 
short-range terms. 



VI. CONCLUSIONS 

This paper has presented a numerical implementation 
of the first-principles envelope-function theory of Ref. 1 
in a model system based on superlattice LDA calcula- 
tions with norm-conserving pseudopotentials. The elec- 
tron density and potential energy of the superlattice were 
approximated by retaining only the linear and quadratic 
response to the heterostructure perturbation. This ap- 
proximation worked well, with a net error of about 2 
meV in GaAs/AlAs and 5 meV in Ino.53Gao.47As/InP. 2 
The principal effect of this error was simply a constant 
shift of the superlattice energy eigenvalues. 

The density and short-range potentials were then ap- 
proximated further using truncated multipole expansions 
(i.e., power series in fc), retaining terms of order k 2 in the 
linear potential and k° in the quadratic potential. This 
had no effect on the macroscopic density and potential 
in bulk, but it generated some additional error (due pri- 
marily to the truncation of the linear density response) 
in a narrow region near the interfaces^ This error was 
confirmed to be negligible for slowly varying envelope 
functions. 

The approximate Hamiltonian was transformed from 
the original plane-wave basis to a Luttinger-Kohn basis 
using zone-center Bloch functions of the reference crys- 
tal. A Luttinger-Kohn unitary transformation was then 
used to eliminate the k • p and potential-energy coupling 
between the A states of interest and the remote B states. 
The resulting basis is material-dependent (due to the 
potential-energy terms) and approximates the position 
dependence of the quasi-Bloch functions in the hetero- 
structure. The perturbation theory of Ref. [l] was ex- 
tended to account for quadratic renormalization of the 



mass and momentum parameters. 

A 7-state {Tisv, T^, ri5 C } envelope-function model 
with linear momentum and mass renormalization 
was shown to give a very good description of the 
Ti5 valence subband structure of GaAs/AlAs and 
Ino.53Gao.47As/InP (001) superlattices, although the 
good results were partly due to a fortuitous cancella- 
tion of errors. Calculations reported elsewhere^ show 
that similar results for Tisv can be obtained from a sim- 
pler 4-state {ri5v,ric} model. A 3-state T 15v model 
gave fairly good results over a more limited energy range 
(although it probably would not work as well in real 
I n o.53Gao.47As/InP superlattices, since the energy gap 
of Ino.53Gao.47As in the model system was significantly 
greater than the experimental value). The primary limi- 
tation of this single-band model is a conflict between the 
need for 0(fc 4 ) bulk terms in order to achieve better ac- 
curacy in the excited states and the sometimes rather se- 
vere plane-wave cutoff needed to avoid spurious solutions 
generated by the 0(fc 4 ) terms. The 3-state model did, 
however, give excellent results for GaAs/Alo^Gao.sAs, 
where the band offset is small enough to satisfy Kohn's 
definition (< 0.1 eV)^ of a shallow perturbation. 

Dipole terms in the quadratic response were found 
to produce interface asymmetry and macroscopic elec- 
tric fields in the no-common-atom Ino.53Gao.47As/InP 
system.- These terms, which have Cn v symmetry, pro- 
duce a significant fraction of the splitting of the quaside- 
generate ground state in such systems. Fitting this split- 
ting to only short-range interface XY mixing terms may 
therefore overestimate the short-range terms and omit 
important physics. 

Numerical results for Ino.53Gao.47As/InP and 
GaAs/AlAs indicate that the linear valence-band 
Rashba coupling parameter is well approximated by 
the bulk effective Lande factor K for cationic pertur- 
bations, but that there is a wide disparity for anionic 
perturbations. Therefore, using bulk magnetoabsorption 
measurements to evaluate interface parameters such 
as the Rashba coupling 5 -^ cannot generally be relied 
upon to provide anything better than a rough order-of- 
magnitude estimate. Of course, the particular numbers 
generated by the present model would likely change 
significantly in a more realistic quasiparticle calculation, 
but the discrepancy between the Rashba and Lande 
coefficients is unlikely to vanish. 

The operator ordering of the effective-mass terms at a 
heterojunction was found to be more complicated than 
in many previous models. Instead of having a single 
von Roos kinetic-energy operator of the form shown in 
Eq. ((T|), perturbation theory yields a linear combina- 
tion of terms with all possible operator orderings. Cer- 
tain terms are larger than others, however. As shown 
by Leibler ; 16 ' 17 to linear order only the BenDaniel- 
Duke operator— Tbd = ^pm~ 1 p and the Gora-Williams 
operator— Tqw = \{ m ~ 1 P 2 + p 2 m~ 1 ) arise. In a simple 
model where the matrix E^ n , in Eq. §§§ is assumed diag- 
onal, the former arises in third-order perturbation theory 
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from the position-dependent energies of remote bands in 
set B, whereas the latter comes from the position depen- 
dence of the bands in set A (see Appendix iDjl. The Zhu- 
Kroemer operator— Tzk = \m~ 1 l 2 p 2 mr 1 / 2 appears as 
one of several terms in quadratic renormalization, and 
the most general von Roos operator does not occur until 
cubic order. 

Actually, with repeated use of the commutator 
[p, f(z)] — —i(df/dz), one can move the momentum oper- 
ators into any desired position. For example, Morrow and 
Brownstein have shown that, upon replacing a — > a + e 
and 7 — * a — e in Eq. |T]), the von Roos operator can be 
rewritten as2i 

1 MlnK)] \ 2 
TvR - Tii -2^{-^z—) ' (22) 

where Th = \m a pmPpm a is the Harrison operator— and 
the second term has the form of a potential energy. But 
one can continue this process indefinitely, for example by 
writing 

T _ d ( 1 d[ln(m")] \ | 1 U[ln(m<*)} \ 2 

dz \2m dz J 2m \ dz J 

Hence, the operator ordering in the effective kinetic en- 
ergy is nothing but an arbitrary convention, as long as 
one takes care to retain all of the effective potential- 
energy terms generated by changing conventions^ 

The effective kinetic-energy operator given by the per- 
turbation theory in this paper does have the advantage 
that the position-dependent functions 9 a (x) appearing in 
it are smooth step-like functions (although it should be 
noted that the position of the step is different for cations 
and anions). One could reduce it to the conventional 
BenDaniel-Duke form, or any other desired form, if one 
were willing to deal with functions having a more com- 
plicated position dependence near the interface. 

This suggests that it may be possible — at least at 
the phenomenological level — to extend the perturbative 
scheme described here to arbitrarily high order in 6^ 
From this perspective, as long as (i) the bulk materi- 
als of the heterostructure are accurately described by an 
effective-mass equation and (ii) the heterostructure per- 
turbation series eventually converges at some finite or- 
der, one can rearrange the operators into some standard 
order, yielding a standardized Hamiltonian with param- 
eters that have a complicated but in principle calculable 
position dependence. 

Of course, this is unlikely to provide a useful first- 
principles calculation method unless the series converges 
at a fairly low order. Fortunately, the examples given 
here demonstrate that linear (in the multiband case) or 
quadratic (in the single-band case) renormalization of the 
momentum and mass parameters is sufficient to achieve 
good results in several typical III-V heterostructures. 
Given that linear-response theory^ has produced accu- 
rate predictions of valence-band offsets in many other 
lattice-matched and lattice-mismatched systems (see the 



bibliography of Ref. [l]), it is likely that the present 
envelope-function theory can be applied successfully in 
many systems too. 
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APPENDIX A: NONLOCAL 
PSEUDOPOTENTIAL 

The nonlocal part of the pseudopotential was handled 
by polynomial fitting in k space. For the bulk refer- 
ence crystal, the entire plane-wave Hamiltonian matrix 
i?(k + G, k + G') = HGG'(k) was evaluated at 57 points 
near k = 0, including k = itself and points of the 
form (100), (110), (111), (200), and (210) x A, where A 
is some specified interval (usually set equal to half the 
magnitude of the smallest superlattice reciprocal-lattice 
vector). It is important to choose a set of fitting points 
with Oh symmetry so that the fitted coefficients main- 
tain the Td symmetry and time-reversal symmetry of the 
Hamiltonian. The Hamiltonian was fitted to a general 
polynomial of order fc 4 with 35 independent coefficients. 
(For general G and G', there is no special symmetry that 
can be used to reduce the number of coefficients.) 

For the linear response, the nonlocal pseudopotential 
T^ nl (k + g, k + g') may be written as V GG , (fc x , k y , k z + 
^(A<? z -(- Ag' z ), Ag z — Ag' z ), in which g is a reciprocal- 
lattice vector of the superlattice and Ag = g — G is 
an integer multiple of (2ir/Nd)z, where N is the num- 
ber of monolayers in the superlattice. The function 
V GG ,(ki, &2, fe, fc/i) was fitted to a general polynomial of 
order k 2 with 15 independent coefficients using a set of 33 
points arranged on a cubic grid. A larger grid was tested 
using fitting polynomials of order k 4 , but the difference 
was negligible so only the simpler method was used. The 
nonlocal pseudopotential is purely linear, so no fitting of 
the quadratic response was necessary. 

APPENDIX B: PERTURBATION THEORY 

This appendix defines a set of functions that offer a 
convenient way to describe operator ordering in fourth- 
order Luttinger-Kohn perturbation theory. These func- 
tions are merely an alternative way of writing the expres- 
sions given on p. 205 of Winkler's monograph.— 

In Luttinger-Kohn perturbation theory,— the total 
Hamiltonian of the system is written as H = Hq + H', 
where Hq has matrix elements (Ho) mm > — E m 5 mm >. The 
states of the unperturbed Hamiltonian Hq are divided 
into a set A containing the states of interest, and a set B 
containing all other states. It is assumed that the ener- 
gies of A and B do not overlap. A unitary transformation 
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H = e~ s He s is used to eliminate the coupling between 
A and B to any desired order in the perturbation H' . 

The notation used here is defined as follows. Mab is 
the block of the matrix M that has rows in set A and 
columns in set B. The matrix G is defined by 

(GUb)„„' = [E n - En')^ 1 , Gba = (Gab) t , (Bl) 

where T denotes the transpose. A dot indicates element- 
by-element multiplication of congruent matrices: 

(A ■ B) nn i = A nn i B nn i , (B2) 

whereas juxtaposition denotes ordinary matrix multipli- 
cation. 

The renormalized Hamiltonian H for states in A is 
given to fourth order in H' by-2& 

Haa = H A A + Pi{H', H') + P 3 (H', H\ H 1 ) 

+ P i (H',H',H',H r ), (B3) 

where the functions P 2 , P3, and P4 are defined by 

P 2 (H\H 2 ) = \[{H\ B ■ Gab)H 2 B a + H\ b {G B a ■ H 2 BA )}, 

(B4) 

P 3 (H\H 2 ,H 3 ) = \{[{H\ B ■ G AB )H 2 BB ] ■ G AB }H BA 
+ \ h ab{ g ba ■ [H BB {G BA ■ H BA )}} 
- j{[ h aa{ g ab ■ H AB )] ■ Gab}H ba 

-\Hab{Gba-[(H 2 ba -G B a)H aa ]}, (B5) 

PiiH 1 , H 2 , H 3 , H 4 ) = 

= \ [( h aa{ g ab ■ [H A a( g ab ■ H AB )]}) ■ Gab\H ba 
+ \ h ab[ g ba ■ ({[{H BA ■ G B a)H aa \ ■ Gba}H aa )] 
~ \ h ab[ g ba ■ {{[H bb (G B a ■ H'ba)] ■ g ba}H aa )] 

- \ h ab[ g ba ■ (H bb {G B a ■ [{H BA ■ g ba)HaaW)} 

- \ [{ h aa{ g ab ■ [{H A b ■ g ab)H bb ]}) ■ g ab]H ba 

- \ [({[ h aa( g ab ■ H AB )] ■ Gab}H bb ) ■ Gab]H ba 

- \ {{[( h ab ■ g ab)H ba \( g ab ■ H AB )} ■ g ab)H ba 

- \ h ab{ g ba ■ {{H B a ■ g ba)[H A b( g ba ■ H BA )]}) 

- \{[{ h ab ■ g ab){H ba ■ g ba)H A b\ ■ g ab}H ba 

- \{\ h ab{ g ba ■ H ba ){Gab ■ H AB )] ■ G A b}H ba 

- \H ab {G B a ■ [H ba ( g ab ■ H ab ){Gba ■ H BA )}} 

- \ h ab{ g ba ■ [( h ba ■ g ba)(H ab ■ g ab)H ba }} 
+ ji^ab • g ab){H ba ■ g ba){H A b ' g ab)H ba 
+ J4 H ab( g ba ■ H ba ){Gab ■ H ab ){Gba ■ H BA ) 
+ \ ( h ab • g ab)(H ba ■ Gba)H ab (G B a ■ H BA ) 
+ \( h ab • g ab)H ba ( g ab ■ H AB )(G BA ■ H BA ) 

+ h[({K H AB ■ g ab)H bb ] ■ Gab}H bb ) ■ G A b\H ba 

+ \h\ b [Gba ■ (H BB {G BA ■ [H BB {G BA ■ H BA )]})]. 

(B6) 



This way of expressing Haa is particularly useful when 
H ' is a sum of operators that do not commute, and one 
wishes to keep track of the order of the various terms. In 
the present case, H' is a sum of k • p terms and potential- 
energy matrix elements. Examples are given in the ap- 
pendices below. 

APPENDIX C: BULK RENORMALIZATION 

In terms of the functions defined in Appendix [B] the 
renormalized coefficients of order k 2 , fc 3 , and k A in the 
bulk reference Hamiltonian for set A are given by-i 

D i L = & i A + P 2 (7r i ,7r j ), (CI) 

G Ta = G fA + P 2^\^ k ) + P 2 {b^^ k ) + P^ l 1 ^^ k ), 

(C2) 

+ P 3 (D ij ,ir k ,ir l )+P 3 (ir l , D ik ,n l ) + P 3 (ir l , n j ,D kl ) 

+ P 2 (D l \D kl ) +P 2 (7T\& kl ) +P 2 (C l ' k ,TT l ). (C3) 

Here a tilde denotes a quantity before renormalization, 
which is obtained by fitting the reference Hamiltonian to 
a polynomial of order fc 4 (as described in Appendix [A| . 
The tilde is omitted on ir l because it does not change 
under renormalization. 

APPENDIX D: LINEAR RENORMALIZATION 

The terms in the renormalized A Hamiltonian that are 
linear in 9 a are given by^ 

*AA = *AA + P2(«\E a ), (Dl) 

7r%A = «AA + P2(E a ,n>), (D2) 

d 7a = Da J a + P ^ a \ ^) + P 2 (E a ,D^) 

+ P 3 (E a ,7T\^), (D3) 

D AA = Daa + P2(K la ,* 3 ) + P 2 {*\* a3 ) 

+ P 3 (7T\E a ,^), (D4) 

d Ta = DTa + p 2 . *** ) + P2 13 , E a ) 

+ P 3 (ir\Tr\E a ). (D5) 

These are the same as the expressions given in Appendix 
C of Ref. G], although written in a different notation. 
Once again, a tilde denotes a quantity before renormal- 
ization, which is obtained from a multipole expansion of 
the linear density and short-range potentials (Sec. Inland 
Ref. 0) and from fitting the linear nonlocal potential to 
a polynomial of order k 2 (Appendix [A"]) . The tilde is 
omitted on E a because it does not change under renor- 
malization. 
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APPENDIX E: QUADRATIC 
RENORMALIZATION 

Perturbative renormalization also generates terms that 
are quadratic in 9 a . The only term included in the Hamil- 
tonian of Ref. Q] was the renormalized potential 

Ef A = P 2 (E a ,E0). (El) 

Some of the present calculations also include quadratic 
renormalization of the momentum matrix: 

(E2) 

<A = P2{E a ^) + P 2 (vP\ EP) + P 3 (E a ,7r\E^), 

(E3) 

*AA = P * ,E fJ ) + P 2 (it 1 ,E a ?) + P 3 (w i ,E a ,EP), 

(E4) 

and of the effective masses: 

D°$ = P 2 (E a ,D^) + P 2 (E a ^m) 
+ P 3 (E a ,E^D^) + P 3 (E a ,w 0i ,^) + P 3 (E^ , tt* , n j ) 

+ P i (E a ,E^*\^), (E5) 

D AA = P 2 , ^ ) + P 3 (^ a ,E^ ) 

+ P 3 (ir\E a ,n^)+P 3 (Tr\E a P,n3) 

+ P 4 (w\E a ,E^wi), (E6) 

D AA = P2{D l3a ,E p ) + P 2 (D l \E afi ) 
+ P 3 (D^,E a ,E^) + P 3 (tt\ n ja ,EP) + P 3 (n\TT^E a ' 3 ) 

+ P i (n\^,E a ,E^), (E7) 

= P 2 (Tr ai ,n i) + P 2 (E a ,D^) 
+ P 3 (E a 1 Ti\^)+P 3 {E a 1 tt* 5 , ir j ) + P 3 (# ai , E? , tp? ) 

+ P i (E a ,w i ,E ,7r j ), (E8) 

Djff = P 2 (^ a ,n^) + P 2 {b ia \Ef>) 
+ P 3 (w ia y,EP) + P 3 (n l , n aj ,E?) + P 3 (n l , E a , n^) 

+ P i {n\E a ,TT\E l3 ) : (E9) 

If*W = P 2 (w ai ,w^) + P 2 (D ai i,EP) + P 2 (E a ,D 1 ^) 
+ P 3 (n m ,^,E^) + P 3 (E a ,D^,EP) + P 3 (E a ,iT\ n^) 

+ P 4 (E a ,ir\Tr j ,E 13 ). (E10) 

In these expressions, several approximations are used. No 
fitting of the momentum and mass terms in the origi- 
nal quadratic Hamiltonian was performed here; conse- 
quently, the unrenormalized parts are set to zero. Also, 



the term E a " is the mean unrenormalized diatomic short- 
range potential summed over all diatomic perturbations 
with the given values of a and f3. 

Some of the calculations also include renormalized 
short-range potential terms of order 9 3 and (9 4 , which 
were approximated using the expressions 

E AA = P2(E ap ,E^) + P 2 (E a , E^) + P 3 (E a , E?, ET), 

(Ell) 

E AA S = p 2{E aP , E^ 5 ) + P a (E a ,E , EP S ) 

+ P 3 (E a ,E^,E s ) + P 3 (E a/3 , E~f , E 5 ) 

+ P i (E ot ,E l *,EP,E 5 ). (E12) 

Corrections of the same order arising from the long-range 
terms in the Hartree potential were neglected. For all of 
the numerical examples treated here, the 0(9 A ) terms in 
Eq. (|E12[) were found to be negligible. 

It should be emphasized that the results presented here 
do not provide a fully consistent perturbation scheme 
according to the criteria given by Takhtamirov and 
Volkov^ in which the mean kinetic energy of the states 
of interest is assumed to be comparable to the hetero- 
structure potential-energy perturbation. According to 
this scheme, if one includes the 0(k 2 9 2 ) and O(k 9 3 ) 
terms shown here, one should also include terms of order 
k 4 9 1 and k 6 9°. 

However, since these require the use of sixth-order 
perturbation theory, such terms were judged to be not 
worth the effort in a preliminary investigation of this na- 
ture. Therefore, the results obtained by adding only the 
0(k 2 9 2 ) and O(k 9 3 ) terms are not expected to be valid 
for kinetic energies covering the full range of the band 
offset, but only for kinetic energies small in comparison 
to the band offset. This is indeed what was found in the 
numerical calculations of Sec. [Vj Likewise, the O(k A 9 ) 
terms were found to be less important for states of small 
kinetic energy. These results are merely a reflection of 
the fact that a wide quantum well, unlike a hydrogenic 
impurity^ does have states in which the mean kinetic en- 
ergy is small compared to the mean perturbing potential. 
Thus, in this sense the theory of low-energy excitations 
of wide quantum wells is actually simpler than the corre- 
sponding theory of shallow impurities, because terms of 
high order in k are of lesser importance. 

APPENDIX F: QUADRATIC RESPONSE 

The method used here to handle the quadratic re- 
sponse differs slightly from that of Ref. [l|. The quadratic 
potential response is given by Eq. (3.14) of Ref. Q]: 

^ 2 >(x,x') = £'*a*R' A^,(x,x'), (Fl) 

q,R a',R' 

which has the same form as the quadratic density re- 
sponse in Eq. (|3c)l . The translation symmetry of the 
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reference crystal allows A/ug^, (x, x') to be written as 

A«g£,(x,x / ) = A« aa '' R '- R (x-R QQS x'-R QQ 0, (F2) 
in which R QCt < is the midpoint of the two atoms: 

Rem' = 7(R + R' + T Q + T Q ' ) . (F3) 

In Eq. (|F2[) , the coordinate reference is taken to be R aa ', 
whereas in Eq. (3.17) of Ref. [l|, it was chosen to be R a = 
R + t q . The Fourier transform of Eq. (|F1[) is 1 

y (2) (k,k') = JV ^ 6» aa ' R '(k-k')Ai; Qa ' R '(k,k'), 

a,a',R' 

(F4) 

where 

r' R '(k) = ^«' +R ' e " MR ° +(R ' +v "° )/2 ' 

(F5) 



The coordinate reference (|F3|1 is arbitrary, but it is some- 
times more convenient for analyzing symmetry properties 
than the choice used in Ref. [fl. 



R 



In the LDA model used here, the quadratic potential 
is local [i.e., Av aa ' R ' (k, k') = Aw QQ ' R '(k - k')] since 
the ionic pseudopotential is purely linear. The diatomic 
potential Av aa R (q) was approximated for small q by 
keeping only the dipole and quadrupole terms in the elec- 
tron density and the I = term in the power-series ex- 
pansion (TjQ) of the short-range potential (which in the 
quadratic case consists only of the exchange-correlation 
potential). For q near G ^ 0, only the I = terms were 
retained in both the density and short-range potential; 
see Sec. VI B of the preceding paper— for further details. 
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